For the potential U we may simply write a sum of terms of the form - Gm_i / r_i, where G is the gravitational constant, m_i the masses of the objects and r_i the distances to the field point. The equation then yields the gravitational field strength. We place the origin at the center of the Earth, and the Moon at a distance d on the y-axis.

Notice that the expression for U in the definitions segment contains square brackets as well as parentheses. Although these brackets have the same effect as the common ones, clever use of pairs of them may improve the readability of an expression.

The contour plot of U yields curves which all lie close to the Earth. In order to study the function over the full domain we request a logarithmic plot of -U, rather than the direct value.

The logarithmic plot of the magnitude of g (Figure 2) demonstrates how the Moon perturbs the Earth's gravitational field. In fact, there is a point close to the Moon where the field even becomes zero.

This is also evident from the vector plot (Figure 3) , where the field direction at that point seems to be undetermined.

Since a field normally has a direction at each point in space, we may also request a contour plot of that direction. This we do in the next plot (not shown), just inverting the equation and then converting from radians to degrees. This leaves an ambiguity concerning the sign of the direction, but in this case this does not cause serious confusion. We can do better with the functions available, but extracting the sign is slightly more complicated.

Planets of Finite Size

In the next descriptor we zoom the fields around each planet. In order to sketch the planets themselves we cut out the areas occupied by the Earth and the Moon. This is also a better presentation in the sense that our potential strictly applies only to the region outside the massive objects.

In order to describe a boundary we have encircled the outside in a positive sense (counter-clockwise). In other words, in walking around the boundary we have previously kept the domain on our left. Cutting out a portion of the domain we still trace the boundary with the domain on our left, which means that we walk in a negative sense. For drawing circles we use the command arc, and the usage should be obvious from the file text.

The grid plot below shows that the program divides the area into smaller triangles as you approach a relatively small object, in this case the planets. This is useful since we are going to inspect the behavior of the gravitational field at close range.

The angular perturbation can be seen even more clearly from the elevation plot (not shown here), which displays the angle of the field as a function of the radial distance from the Moon. Evidently the distortion is only about 0.1 degrees at the surface of the Moon. Thus, an astronaut should not find it difficult to remain on his feet.

Figure 7 is a zoomed (and rotated) surface plot of the magnitude of the gravitational field close to the Moon. We have used a logarithmic scale for the field strength in order to see the Moon's near field as well as the sharp minimum (in fact a zero point) at some distance toward the Earth.

Divergence of g

Before leaving the gravitational field we should investigate whether the divergence is equal to zero. For this purpose we add the following definitions to a copy of grav2.pde and make the corresponding contour plots.

laplace_xy=dxx(U)+dyy(U)

laplace_rz=(1/x)*dx[x*dx(U)]+dyy(U)

To the reader who knows that the divergence of a gravitational field in free space vanishes, the plot of laplace_xy may come as a shock. We are now in three dimensions, however, although the plots are in a plane through the two planetary objects. If we had included the term dzz(U), which is not allowed in the present versions of the software, we would in fact have obtained a vanishingly small value.

An alternative solution would be to change to cylindrical coordinates, ( and z, which we shall have many occasions to use. Plotting laplace_rz we would recognize the irregular pattern of zero contours which indicates a vanishingly small function.

Towards the Center of the Earth

The gravitational field of course does not stop at the surface of the Earth, but penetrates through all of the planet. Newton already showed that the field at any radius below the surface is equivalent to the attraction from all the mass contained within that radius, concentrated at the center. This does not mean that we may compute the potential by summation over internal masses only. What we can do is to derive the potential in this region from Newton's expression for g, assuming the Earth's density to be constant. It is easy to show, by taking the derivative, that the following expression for U yields the correct gravitational field on the inside, and that the potential takes the previous value at the Earth's radius r₁.

We now study this potential using much the same descriptor as before, as suggested below.

Here, we have expressed the radial coordinate rad in terms of Cartesian x and y. The plots show that the force on an object in this region is proportional to the radius and directed towards the center. The surface plot below displays the magnitude of g.

In addition we introduce the Laplacian of U, using x temporarily as the radial coordinate in a polar system. The elevation plot of laplace_U for y=0 demonstrates that the result is a constant value. This illustrates that the potential obeys the partial differential equation for gravitation (Poisson's equation)

and you can easily verify that the value we obtain for the Laplacian is in accord with the uniform density assumed.